R Notebook

Colonrectal cancer project: reproduce GSE149224 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7486130/

count <- read.delim(gzfile("D:\\scRNAseq\\GSE149224_RSH.all.counts.txt.gz"), sep = " ")
count

NA

# count_2 <- read.delim(gzfile("D:\\scRNAseq\\GSE149224_RSH.all.txt.gz"), sep = " ")

annot <- read.delim(gzfile("D:\\scRNAseq\\GSE149224_meta.information.csv.gz"), sep = ",")
annot <- annot %>% unite(batch_gid_dose, c(batch, df.gid,dose), sep = "_", remove = FALSE)
annot

unique(annot$df.gid)

[1] "HCT116" "SW480"  "RKO"   

gene_names <- rownames(count)
unique(length(gene_names))

[1] 23768

Convert gene symbols to entrez

library("EnsDb.Hsapiens.v86")

Loading required package: ensembldb
Loading required package: GenomicFeatures
Loading required package: AnnotationFilter

Attaching package: 'ensembldb'

The following object is masked from 'package:dplyr':

    filter

The following object is masked from 'package:stats':

    filter

mapIds_entrezid <- mapIds(EnsDb.Hsapiens.v86, keys = gene_names, keytype = "GENENAME", columns = c("ENTREZID"))

Warning: Unable to map 1936 of 23768 requested IDs.

mapIds_GENEID <- mapIds(EnsDb.Hsapiens.v86, keys = gene_names, keytype = "GENENAME", columns = c("GENEID"))

Warning: Unable to map 1936 of 23768 requested IDs.

crc_dat <- SingleCellExperiment(
  assays = list(counts = as.matrix(count)),
  colData = annot
)
crc_dat

class: SingleCellExperiment 
dim: 23768 11126 
metadata(0):
assays(1): counts
rownames(23768): RP11-34P13.7 AL627309.1 ... KDM5D EIF1AY
rowData names(0):
colnames(11126): X112385.TCCGTCTTGATC X112385.AGTTTCCCTGGN ... X120386.CGATGGCCGAAT
  X120386.AACGGTGCATAA
colData names(5): X batch_gid_dose batch df.gid dose
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

After we create our single cell object, we remove variables that we don’t use later for memnory saving purpose since rnaseq file are huge in most cases. We use rm() command to do that job

rm(c("count", "annot"))

Error in rm(c("count", "annot")) : 
  ... must contain names or character strings

We, we create a log2 transform assay for our scRNA object

assay(crc_dat, "log2counts") <- log2(counts(crc_dat) + 1)

We take a peek on the coldata information. coldata has information of the RNA transcript and related information to the experiment of samples

crc_dat@colData

DataFrame with 11126 rows and 5 columns
                                       X   batch_gid_dose     batch      df.gid      dose
                             <character>      <character> <integer> <character> <integer>
X112385.TCCGTCTTGATC 112385-TCCGTCTTGATC  112385_HCT116_0    112385      HCT116         0
X112385.AGTTTCCCTGGN 112385-AGTTTCCCTGGN   112385_SW480_0    112385       SW480         0
X112385.TGTTCCCATTCN 112385-TGTTCCCATTCN   112385_SW480_0    112385       SW480         0
X112385.CGGGACACTCCC 112385-CGGGACACTCCC  112385_HCT116_0    112385      HCT116         0
X112385.ACCGGGTTACGG 112385-ACCGGGTTACGG  112385_HCT116_0    112385      HCT116         0
...                                  ...              ...       ...         ...       ...
X120386.TACAAGATATAC 120386-TACAAGATATAC 120386_SW480_200    120386       SW480       200
X120386.CGCGATCAAATC 120386-CGCGATCAAATC 120386_SW480_200    120386       SW480       200
X120386.AGACTGGCTTTT 120386-AGACTGGCTTTT 120386_SW480_200    120386       SW480       200
X120386.CGATGGCCGAAT 120386-CGATGGCCGAAT 120386_SW480_200    120386       SW480       200
X120386.AACGGTGCATAA 120386-AACGGTGCATAA 120386_SW480_200    120386       SW480       200

Next, we look at the mean and total counts for each gene by creating those columns in colData of our single cell object

colData(crc_dat)$mean <- colMeans(counts(crc_dat))
colData(crc_dat)$total_count <- colCounts(counts(crc_dat))

# Making another assay cpm
assay(crc_dat, "cpm") <- counts(crc_dat)/crc_dat$total_count/1e6 # Count for cell per million

We transform column that need to be factor data into factor

colData(crc_dat)

DataFrame with 11126 rows and 7 columns
                                       X   batch_gid_dose    batch   df.gid      dose      mean total_count
                             <character>         <factor> <factor> <factor> <integer> <numeric>   <integer>
X112385.TCCGTCTTGATC 112385-TCCGTCTTGATC  112385_HCT116_0   112385   HCT116         0  0.802255        1815
X112385.AGTTTCCCTGGN 112385-AGTTTCCCTGGN  112385_SW480_0    112385   SW480          0  0.180032        1265
X112385.TGTTCCCATTCN 112385-TGTTCCCATTCN  112385_SW480_0    112385   SW480          0  0.365533        1549
X112385.CGGGACACTCCC 112385-CGGGACACTCCC  112385_HCT116_0   112385   HCT116         0  0.140483        1001
X112385.ACCGGGTTACGG 112385-ACCGGGTTACGG  112385_HCT116_0   112385   HCT116         0  0.117721         911
...                                  ...              ...      ...      ...       ...       ...         ...
X120386.TACAAGATATAC 120386-TACAAGATATAC 120386_SW480_200   120386    SW480       200 0.0377819         428
X120386.CGCGATCAAATC 120386-CGCGATCAAATC 120386_SW480_200   120386    SW480       200 0.0457758         501
X120386.AGACTGGCTTTT 120386-AGACTGGCTTTT 120386_SW480_200   120386    SW480       200 0.0546112         135
X120386.CGATGGCCGAAT 120386-CGATGGCCGAAT 120386_SW480_200   120386    SW480       200 0.0359307         518
X120386.AACGGTGCATAA 120386-AACGGTGCATAA 120386_SW480_200   120386    SW480       200 0.0359727         462

Now, we will take a look at the violin plot of total count per batch and dose and df.gid

ggplot(data = cell_info, aes(x = batch_gid_dose, y = total_count)) +
  geom_violin(fill = 'brown') + theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

Normally, if your gene names are in ENS format, we may need to convert to symbols so we can easily to do gene annotation later for gene set analysis or pathway analysis. However, we already provided genes in symbol format, we don’t need to do further

Let’s check if we can find mitochondrial proteins in the newly annotated symbols. Mitochrondial is an essential part of QC since if its proportion is too large (more than 10%, usually) then the data is not good to use

Our data given symbols for gene names, we may want to convert to other types for later use if needed. Here, I converted to entrezid and geneid

rowData(crc_dat)$SYMBOL <- gene_names
rowData(crc_dat)$entrezid <- mapIds_entrezid
rowData(crc_dat)$geneid <- mapIds_GENEID

We can check mitochrondial in our sequence

mito_names <- grep("^MT-",rowData(crc_dat)$SYMBOL,value = T)
mito_names

 [1] "MT-TF"   "MT-RNR1" "MT-TV"   "MT-RNR2" "MT-TL1"  "MT-ND1"  "MT-TI"   "MT-TQ"   "MT-TM"   "MT-ND2"  "MT-TW"  
[12] "MT-TA"   "MT-TN"   "MT-TC"   "MT-TY"   "MT-CO1"  "MT-TS1"  "MT-TD"   "MT-CO2"  "MT-TK"   "MT-ATP8" "MT-ATP6"
[23] "MT-CO3"  "MT-TG"   "MT-ND3"  "MT-TR"   "MT-ND4L" "MT-ND4"  "MT-TH"   "MT-TS2"  "MT-TL2"  "MT-ND5"  "MT-ND6" 
[34] "MT-TE"   "MT-CYB"  "MT-TT"   "MT-TP"  

And find ribosomal proteins if you want (high proportion of ribosomal proteins may indicate that cells are dying/dead)

• Similar command to find ribosomal proteins (which start with RPL or RPS) works as expected:

ribo_names <- grep("^RP[LS]",rowData(crc_dat)$SYMBOL,value = T)
ribo_names

  [1] "RPL22"          "RPL11"          "RPS6KA1"        "RPS8"          
  [5] "RPL5"           "RPS27"          "RPS10P7"        "RPS6KC1"       
  [9] "RPS24"          "RPLP2"          "RPL27A"         "RPS13"         
 [13] "RPS6KA4"        "RPS6KB2"        "RPS3"           "RPS25"         
 [17] "RPS26"          "RPL41"          "RPL6"           "RPLP0"         
 [21] "RPL21"          "RPL10L"         "RPS29"          "RPL36AL"       
 [25] "RPS6KL1"        "RPS6KA5"        "RPS27L"         "RPL4"          
 [29] "RPLP1"          "RPS17"          "RPS17L"         "RPL3L"         
 [33] "RPS2"           "RPS15A"         "RPL13"          "RPL26"         
 [37] "RPL23A"         "RPL23"          "RPL19"          "RPL27"         
 [41] "RPS6KB1"        "RPL38"          "RPL17-C18orf32" "RPL17"         
 [45] "RPS15"          "RPL36"          "RPS28"          "RPL18A"        
 [49] "RPSAP58"        "RPS16"          "RPS19"          "RPL18"         
 [53] "RPL13A"         "RPS11"          "RPS9"           "RPL28"         
 [57] "RPS5"           "RPS7"           "RPS27A"         "RPL31"         
 [61] "RPL37A"         "RPS21"          "RPL3"           "RPS19BP1"      
 [65] "RPL32"          "RPL15"          "RPSA"           "RPL14"         
 [69] "RPL29"          "RPL24"          "RPL22L1"        "RPL39L"        
 [73] "RPL35A"         "RPL9"           "RPL34-AS1"      "RPL34"         
 [77] "RPS3A"          "RPL37"          "RPS23"          "RPS14"         
 [81] "RPL26L1"        "RPS18"          "RPS10-NUDT3"    "RPS10"         
 [85] "RPL10A"         "RPL7L1"         "RPS12"          "RPS6KA2"       
 [89] "RPS20"          "RPL7"           "RPL30"          "RPL8"          
 [93] "RPS6"           "RPL35"          "RPL12"          "RPL7A"         
 [97] "RPS6KA3"        "RPS4X"          "RPS6KA6"        "RPL36A"        
[101] "RPL39"          "RPL10"          "RPS4Y1"        

Quick search for mitochondrial protein ATP8, which is also called MT-ATP8

grep("*ATP8",rowData(crc_dat)$SYMBOL,value = T)

[1] "ATP8B2"  "ATP8A2"  "ATP8B4"  "ATP8B1"  "ATP8B3"  "ATP8A1"  "MT-ATP8"

It is time for us to take a look at more details of the data and do some QC for the dataset.

We do some visualization regarding number of cells, percentage, …

+The following scater functions allow us to add per-cell and per-gene metrics useful for dataset evaluation. Most popular metrics per cell are total number of counts (UMIs), total number of detected genes, total number of mitochondrial counts, percent of mitochondrial counts, etc.

First, we look at the counts and percentage of each groups correspoding to dose of treatments

df_plot <- data.frame(dose = factor(crc_dat@colData$dose), group = factor(crc_dat@colData$batch_gid_dose))

df_plot %>% 
    ggplot(aes(x=dose, fill=group)) + 
    geom_bar() +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    theme(plot.title = element_text(hjust=0.5, face="bold")) +
    ggtitle("NCells")

And the density distribution of each batch for dose, gid and batch

# Visualize the number UMIs/transcripts per cell
data.frame(batch_gid_dose = crc_dat@colData$batch_gid_dose, total_count = crc_dat@colData$total_count) %>%
    ggplot(aes(color=batch_gid_dose, x=total_count, fill= batch_gid_dose)) + 
    geom_density(alpha = 0.2) + 
    scale_x_log10() + 
    theme_classic() +
    ylab("Cell density") +
    geom_vline(xintercept = 500)

Now, we just take a closer look at the percentage of Mitochrondial cells and ribosomal proteins percentage

crc_dat_Mitocell <- perCellQCMetrics(crc_dat,subsets=list(Mito=mito_names))

crc_dat_ribocell <- perCellQCMetrics(crc_dat,subsets=list(Ribo=ribo_names))
crc_dat_feature <- perFeatureQCMetrics(crc_dat)

• If you want to modify the colData by adding those information, we can use the following way

crc_dat <- addPerCellQC(crc_dat, subsets=list(Mito=mito_names))
crc_dat <- addPerCellQC(crc_dat, subsets=list(Ribo=ribo_names))

names(colData(crc_dat))

 [1] "X"                     "batch_gid_dose"        "batch"                
 [4] "df.gid"                "dose"                  "mean"                 
 [7] "total_count"           "sum"                   "detected"             
[10] "subsets_Mito_sum"      "subsets_Mito_detected" "subsets_Mito_percent" 
[13] "total"                 "sum"                   "detected"             
[16] "subsets_Ribo_sum"      "subsets_Ribo_detected" "subsets_Ribo_percent" 
[19] "total"                

And the table of colData

Or we can see them separately if we use the separate results

Mitochrondial percentage

crc_dat_Mitocell

DataFrame with 11126 rows and 6 columns
                           sum  detected subsets_Mito_sum subsets_Mito_detected
                     <numeric> <numeric>        <numeric>             <numeric>
X112385.TCCGTCTTGATC     19068      3946              646                    17
X112385.AGTTTCCCTGGN      4279      2014              212                    15
X112385.TGTTCCCATTCN      8688      2835              236                    14
X112385.CGGGACACTCCC      3339      1564              166                    15
X112385.ACCGGGTTACGG      2798      1409              125                    12
...                        ...       ...              ...                   ...
X120386.TACAAGATATAC       898       561               64                    11
X120386.CGCGATCAAATC      1088       671               91                    11
X120386.AGACTGGCTTTT      1298       161              103                     9
X120386.CGATGGCCGAAT       854       639               58                    12
X120386.AACGGTGCATAA       855       592               94                    12
                     subsets_Mito_percent     total
                                <numeric> <numeric>
X112385.TCCGTCTTGATC              3.38787     19068
X112385.AGTTTCCCTGGN              4.95443      4279
X112385.TGTTCCCATTCN              2.71639      8688
X112385.CGGGACACTCCC              4.97155      3339
X112385.ACCGGGTTACGG              4.46748      2798
...                                   ...       ...
X120386.TACAAGATATAC              7.12695       898
X120386.CGCGATCAAATC              8.36397      1088
X120386.AGACTGGCTTTT              7.93529      1298
X120386.CGATGGCCGAAT              6.79157       854
X120386.AACGGTGCATAA             10.99415       855

Ribosomal percentage

crc_dat_ribocell

DataFrame with 11126 rows and 6 columns
                           sum  detected subsets_Ribo_sum subsets_Ribo_detected
                     <numeric> <numeric>        <numeric>             <numeric>
X112385.TCCGTCTTGATC     19068      3946             3509                    86
X112385.AGTTTCCCTGGN      4279      2014              616                    74
X112385.TGTTCCCATTCN      8688      2835             1685                    81
X112385.CGGGACACTCCC      3339      1564              547                    79
X112385.ACCGGGTTACGG      2798      1409              439                    69
...                        ...       ...              ...                   ...
X120386.TACAAGATATAC       898       561              222                    59
X120386.CGCGATCAAATC      1088       671              198                    61
X120386.AGACTGGCTTTT      1298       161               19                    15
X120386.CGATGGCCGAAT       854       639               77                    35
X120386.AACGGTGCATAA       855       592               88                    44
                     subsets_Ribo_percent     total
                                <numeric> <numeric>
X112385.TCCGTCTTGATC              18.4026     19068
X112385.AGTTTCCCTGGN              14.3959      4279
X112385.TGTTCCCATTCN              19.3946      8688
X112385.CGGGACACTCCC              16.3822      3339
X112385.ACCGGGTTACGG              15.6898      2798
...                                   ...       ...
X120386.TACAAGATATAC             24.72160       898
X120386.CGCGATCAAATC             18.19853      1088
X120386.AGACTGGCTTTT              1.46379      1298
X120386.CGATGGCCGAAT              9.01639       854
X120386.AACGGTGCATAA             10.29240       855

Find low-quality cells

We then identify cells that are outliers for the various QC metrics, based on the median absolute deviation (MAD) from the median value of each metric across all cells. Specifically, a value is considered an outlier if it is more than 3 MADs from the median in the “problematic” direction. This is loosely motivated by the fact that such a filter will retain 99% of non-outlier values that follow a normal distribution.

Identify outliers

• Find cell that has outside the 3MAD from median for mitochrondial

mito_sum_low <- isOutlier(crc_dat_Mitocell$sum, log=TRUE, type="lower")

mito_detect_high <- isOutlier(crc_dat_Mitocell$detected, log=TRUE, type="higher")

• Check the finding outliers

attr(mito_low, "thresholds")

   lower   higher 
250.7705      Inf 

attr(mito_high, "thresholds")

  lower  higher 
   0.00 4794.34 

• Find the cutoff for percentage for mito

mito_per <- isOutlier(crc_dat_Mitocell$subsets_Mito_percent, type="higher")
attr(mito_per, "thresholds")

   lower   higher 
    -Inf 12.35193 

• We remove the outliers and summary what we did

# We take all the ones that does not satistfy criterions by removing low size library, highe detected library and high percentage of mitochrondial 

mito_discard <- mito_detect_high | mito_sum_low | mito_per 

# Summarize the number of cells removed for each reason.
DataFrame(LibSize=sum(mito_sum_low), NExprs=sum(mito_detect_high),
     MitoProp=sum(mito_per), Total=sum(mito_discard))

• Another way to do the tasks above

reasons

DataFrame with 11126 rows and 4 columns
          low_lib_size   low_n_features high_subsets_Mito_percent   discard
      <outlier.filter> <outlier.filter>          <outlier.filter> <logical>
1                FALSE            FALSE                     FALSE     FALSE
2                FALSE            FALSE                     FALSE     FALSE
3                FALSE            FALSE                     FALSE     FALSE
4                FALSE            FALSE                     FALSE     FALSE
5                FALSE            FALSE                     FALSE     FALSE
...                ...              ...                       ...       ...
11122            FALSE            FALSE                     FALSE     FALSE
11123            FALSE            FALSE                     FALSE     FALSE
11124            FALSE             TRUE                     FALSE      TRUE
11125            FALSE            FALSE                     FALSE     FALSE
11126            FALSE            FALSE                     FALSE     FALSE

• Do some visualization

Now, we take a look at the cells that we may remove by some features

• By df.gid

 plotColData(crc_dat, x="batch", y="sum", colour_by="discard",
        other_fields="df.gid") + facet_wrap(~df.gid, nrow = 4) + 
        scale_y_log10() + ggtitle("Total count")

• By dose

 plotColData(crc_dat, x="batch", y="sum", colour_by="discard",
        other_fields="dose") + facet_wrap(~dose, nrow = 4) + 
        scale_y_log10() + ggtitle("Total count")

Subset the data base on discard criteria

crc_dat_filter

class: SingleCellExperiment 
dim: 23768 10782 
metadata(0):
assays(3): counts log2counts cpm
rownames(23768): RP11-34P13.7 AL627309.1 ... KDM5D EIF1AY
rowData names(3): SYMBOL entrezid geneid
colnames(10782): X112385.TCCGTCTTGATC X112385.AGTTTCCCTGGN ...
  X120386.CGATGGCCGAAT X120386.AACGGTGCATAA
colData names(20): X batch_gid_dose ... total discard
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

• Before the discard

crc_dat

class: SingleCellExperiment 
dim: 23768 11126 
metadata(0):
assays(3): counts log2counts cpm
rownames(23768): RP11-34P13.7 AL627309.1 ... KDM5D EIF1AY
rowData names(3): SYMBOL entrezid geneid
colnames(11126): X112385.TCCGTCTTGATC X112385.AGTTTCCCTGGN ...
  X120386.CGATGGCCGAAT X120386.AACGGTGCATAA
colData names(20): X batch_gid_dose ... total discard
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

saveRDS(crc_dat, file = "D:\\scRNAseq\\crc_dat.rds")
saveRDS(crc_dat_filter, file = "D:\\scRNAseq\\crc_dat_filter.rds")

Load dataset that we save above

Normalization by using librarySizeFactors from scater package

library(scater)

Loading required package: scuttle

lib.crc_dat_filter <- librarySizeFactors(crc_dat_filter)
summary(lib.crc_dat_filter)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1209  0.4921  0.7765  1.0000  1.2285  8.8894 

Histogram of the normalization

hist(log10(lib.crc_dat_filter), xlab="Log10[Size factor]", col='grey80')

Normalization by deconvolution

library(scran)
set.seed(100)
clust.crc_dat_filter <- quickCluster(crc_dat_filter) 
table(clust.crc_dat_filter)

clust.crc_dat_filter
   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
1369 1034 1077 1044  725  680  887  590  981  606  445  583  133  264  146  218 

deconv.crc_dat_filter <- calculateSumFactors(crc_dat_filter, cluster=clust.crc_dat_filter)
summary(deconv.crc_dat_filter)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.1047  0.4827  0.7611  1.0000  1.2207  8.5547 

We use a pre-clustering step with quickCluster() where cells in each cluster are normalized separately and the size factors are rescaled to be comparable across clusters. This avoids the assumption that most genes are non-DE across the entire population - only a non-DE majority is required between pairs of clusters, which is a weaker assumption for highly heterogeneous populations. By default, quickCluster() will use an approximate algorithm for PCA based on methods from the irlba package. The approximation relies on stochastic initialization so we need to set the random seed (via set.seed()) for reproducibility.

plot(lib.crc_dat_filter, deconv.crc_dat_filter, xlab="Library size factor",
    ylab="Deconvolution size factor", log='xy', pch=16,
    col=as.integer(factor(crc_dat_filter$batch)))
abline(a=0, b=1, col="red")

Next, we look at the logtransform

• We apply the logNormCounts function from scater package

crc_dat_filter <- logNormCounts(crc_dat_filter)
assayNames(crc_dat_filter)

[1] "counts"     "log2counts" "cpm"        "logcounts" 

crc_dat_filter@assays@data$log2counts[1:10,1:10]

               X112385.TCCGTCTTGATC X112385.AGTTTCCCTGGN X112385.TGTTCCCATTCN X112385.CGGGACACTCCC
RP11-34P13.7                      0                    0                    0                    0
AL627309.1                        0                    0                    0                    0
RP11-34P13.13                     0                    0                    0                    0
AP006222.2                        0                    0                    0                    0
RP5-857K21.3                      0                    0                    0                    0
RP11-206L10.3                     0                    0                    0                    0
RP11-206L10.2                     0                    0                    0                    0
RP11-206L10.9                     0                    0                    0                    0
RP11-206L10.10                    0                    0                    0                    0
FAM87B                            0                    0                    0                    0
               X112385.ACCGGGTTACGG X112385.GTCTTTTCAAGC X112385.GAGCAACGTCAN X112385.CCTAATGCGGTN
RP11-34P13.7                      0                    0                    0                    0
AL627309.1                        0                    0                    0                    0
RP11-34P13.13                     0                    0                    0                    0
AP006222.2                        0                    0                    0                    0
RP5-857K21.3                      0                    0                    0                    0
RP11-206L10.3                     0                    0                    0                    0
RP11-206L10.2                     0                    0                    0                    0
RP11-206L10.9                     0                    0                    0                    0
RP11-206L10.10                    0                    0                    0                    0
FAM87B                            0                    0                    0                    0
               X112385.ACCATACACGCN X112385.GCCGTGGGCTGN
RP11-34P13.7                      0                    0
AL627309.1                        0                    1
RP11-34P13.13                     0                    0
AP006222.2                        0                    0
RP5-857K21.3                      0                    0
RP11-206L10.3                     0                    0
RP11-206L10.2                     0                    0
RP11-206L10.9                     0                    0
RP11-206L10.10                    0                    0
FAM87B                            0                    0

crc_dat_filter@assays@data$logcounts[1:10,1:10]

               X112385.TCCGTCTTGATC X112385.AGTTTCCCTGGN X112385.TGTTCCCATTCN X112385.CGGGACACTCCC
RP11-34P13.7                      0                    0                    0                    0
AL627309.1                        0                    0                    0                    0
RP11-34P13.13                     0                    0                    0                    0
AP006222.2                        0                    0                    0                    0
RP5-857K21.3                      0                    0                    0                    0
RP11-206L10.3                     0                    0                    0                    0
RP11-206L10.2                     0                    0                    0                    0
RP11-206L10.9                     0                    0                    0                    0
RP11-206L10.10                    0                    0                    0                    0
FAM87B                            0                    0                    0                    0
               X112385.ACCGGGTTACGG X112385.GTCTTTTCAAGC X112385.GAGCAACGTCAN X112385.CCTAATGCGGTN
RP11-34P13.7                      0                    0                    0                    0
AL627309.1                        0                    0                    0                    0
RP11-34P13.13                     0                    0                    0                    0
AP006222.2                        0                    0                    0                    0
RP5-857K21.3                      0                    0                    0                    0
RP11-206L10.3                     0                    0                    0                    0
RP11-206L10.2                     0                    0                    0                    0
RP11-206L10.9                     0                    0                    0                    0
RP11-206L10.10                    0                    0                    0                    0
FAM87B                            0                    0                    0                    0
               X112385.ACCATACACGCN X112385.GCCGTGGGCTGN
RP11-34P13.7                      0            0.0000000
AL627309.1                        0            0.4583876
RP11-34P13.13                     0            0.0000000
AP006222.2                        0            0.0000000
RP5-857K21.3                      0            0.0000000
RP11-206L10.3                     0            0.0000000
RP11-206L10.2                     0            0.0000000
RP11-206L10.9                     0            0.0000000
RP11-206L10.10                    0            0.0000000
FAM87B                            0            0.0000000

Dimention reduction analysis

• We perform the PCA on the log-normalized expression values using the runPCA() function from scater. By default, runPCA() will compute the first 50 PCs and store them in the reducedDims() of the output SingleCellExperiment object, as shown below. Here, we use only the top 2000 genes with the largest biological components to reduce both computational work and high-dimensional random noise. In particular, while PCA is robust to random noise, an excess of it may cause the earlier PCs to capture noise instead of biological structure (Johnstone and Lu 2009). This effect can be mitigated by restricting the PCA to a subset of HVGs

library(scran)
top.crc <- getTopHVGs(crc_dat_filter, n=2000)

library(scater)
set.seed(100) # See below.
crc_dat_filter <- runPCA(crc_dat_filter, subset_row=top.crc) 
reducedDimNames(crc_dat_filter)

[1] "PCA"

dim(reducedDim(crc_dat_filter, "PCA"))

[1] 10782    50

How many PCs would we take? We can use the elbow method to find the number of PCs.

library(PCAtools)

Loading required package: ggrepel

Attaching package: ‘PCAtools’

The following object is masked from ‘package:scran’:

    parallelPCA

The following objects are masked from ‘package:stats’:

    biplot, screeplot

percent.var <- attr(reducedDim(crc_dat_filter), "percentVar")
plot(percent.var, log="y", xlab="PC", ylab="Variance explained (%)")

chosen.elbow <- PCAtools::findElbowPoint(percent.var)
chosen.elbow

[1] 4

plot(percent.var, xlab="PC", ylab="Variance explained (%)")
abline(v=chosen.elbow, col="red")

Next we will run UMAP on the top of PCs

Cluster by using k-mean method

set.seed(100)
clust.kmeans <- kmeans(reducedDim(crc_dat_filter, "PCA"), centers=15)
table(clust.kmeans$cluster)

   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
 697 1391  709  384 1492 1313  661  219 1067  636  552  498  696  262  205 

library(SingleR)

# celldex has many different type of cell identification, we use DatabaseImmuneCellExpressionData

monaco.ref <- celldex::MonacoImmuneData()

  |                                                                                                               
  |                                                                                                         |   0%
  |                                                                                                               
  |=                                                                                                        |   1%
  |                                                                                                               
  |==                                                                                                       |   2%
  |                                                                                                               
  |===                                                                                                      |   2%
  |                                                                                                               
  |===                                                                                                      |   3%
  |                                                                                                               
  |====                                                                                                     |   4%
  |                                                                                                               
  |=====                                                                                                    |   5%
  |                                                                                                               
  |======                                                                                                   |   5%
  |                                                                                                               
  |======                                                                                                   |   6%
  |                                                                                                               
  |=======                                                                                                  |   7%
  |                                                                                                               
  |========                                                                                                 |   7%
  |                                                                                                               
  |========                                                                                                 |   8%
  |                                                                                                               
  |=========                                                                                                |   9%
  |                                                                                                               
  |==========                                                                                               |   9%
  |                                                                                                               
  |==========                                                                                               |  10%
  |                                                                                                               
  |===========                                                                                              |  10%
  |                                                                                                               
  |===========                                                                                              |  11%
  |                                                                                                               
  |============                                                                                             |  11%
  |                                                                                                               
  |============                                                                                             |  12%
  |                                                                                                               
  |=============                                                                                            |  12%
  |                                                                                                               
  |=============                                                                                            |  13%
  |                                                                                                               
  |==============                                                                                           |  13%
  |                                                                                                               
  |==============                                                                                           |  14%
  |                                                                                                               
  |==================                                                                                       |  17%
  |                                                                                                               
  |==================                                                                                       |  18%
  |                                                                                                               
  |===================                                                                                      |  18%
  |                                                                                                               
  |====================                                                                                     |  19%
  |                                                                                                               
  |=====================                                                                                    |  20%
  |                                                                                                               
  |======================                                                                                   |  21%
  |                                                                                                               
  |=======================                                                                                  |  22%
  |                                                                                                               
  |========================                                                                                 |  23%
  |                                                                                                               
  |=========================                                                                                |  24%
  |                                                                                                               
  |==========================                                                                               |  25%
  |                                                                                                               
  |===========================                                                                              |  26%
  |                                                                                                               
  |============================                                                                             |  26%
  |                                                                                                               
  |============================                                                                             |  27%
  |                                                                                                               
  |=============================                                                                            |  27%
  |                                                                                                               
  |=============================                                                                            |  28%
  |                                                                                                               
  |==============================                                                                           |  28%
  |                                                                                                               
  |==============================                                                                           |  29%
  |                                                                                                               
  |===============================                                                                          |  29%
  |                                                                                                               
  |===============================                                                                          |  30%
  |                                                                                                               
  |================================                                                                         |  31%
  |                                                                                                               
  |=================================                                                                        |  31%
  |                                                                                                               
  |=================================                                                                        |  32%
  |                                                                                                               
  |==================================                                                                       |  32%
  |                                                                                                               
  |==================================                                                                       |  33%
  |                                                                                                               
  |===================================                                                                      |  33%
  |                                                                                                               
  |===================================                                                                      |  34%
  |                                                                                                               
  |====================================                                                                     |  34%
  |                                                                                                               
  |====================================                                                                     |  35%
  |                                                                                                               
  |=====================================                                                                    |  35%
  |                                                                                                               
  |======================================                                                                   |  36%
  |                                                                                                               
  |======================================                                                                   |  37%
  |                                                                                                               
  |=======================================                                                                  |  37%
  |                                                                                                               
  |========================================                                                                 |  38%
  |                                                                                                               
  |========================================                                                                 |  39%
  |                                                                                                               
  |=========================================                                                                |  39%
  |                                                                                                               
  |==========================================                                                               |  40%
  |                                                                                                               
  |===========================================                                                              |  41%
  |                                                                                                               
  |============================================                                                             |  42%
  |                                                                                                               
  |=============================================                                                            |  43%
  |                                                                                                               
  |==============================================                                                           |  44%
  |                                                                                                               
  |===============================================                                                          |  45%
  |                                                                                                               
  |================================================                                                         |  45%
  |                                                                                                               
  |================================================                                                         |  46%
  |                                                                                                               
  |=================================================                                                        |  46%
  |                                                                                                               
  |=================================================                                                        |  47%
  |                                                                                                               
  |==================================================                                                       |  47%
  |                                                                                                               
  |==================================================                                                       |  48%
  |                                                                                                               
  |===================================================                                                      |  48%
  |                                                                                                               
  |===================================================                                                      |  49%
  |                                                                                                               
  |====================================================                                                     |  49%
  |                                                                                                               
  |====================================================                                                     |  50%
  |                                                                                                               
  |=====================================================                                                    |  50%
  |                                                                                                               
  |======================================================                                                   |  51%
  |                                                                                                               
  |======================================================                                                   |  52%
  |                                                                                                               
  |=======================================================                                                  |  52%
  |                                                                                                               
  |=======================================================                                                  |  53%
  |                                                                                                               
  |========================================================                                                 |  53%
  |                                                                                                               
  |========================================================                                                 |  54%
  |                                                                                                               
  |=========================================================                                                |  54%
  |                                                                                                               
  |=========================================================                                                |  55%
  |                                                                                                               
  |==========================================================                                               |  55%
  |                                                                                                               
  |===========================================================                                              |  56%
  |                                                                                                               
  |===========================================================                                              |  57%
  |                                                                                                               
  |============================================================                                             |  57%
  |                                                                                                               
  |=============================================================                                            |  58%
  |                                                                                                               
  |==============================================================                                           |  59%
  |                                                                                                               
  |================================================================                                         |  61%
  |                                                                                                               
  |==================================================================                                       |  63%
  |                                                                                                               
  |====================================================================                                     |  65%
  |                                                                                                               
  |=====================================================================                                    |  66%
  |                                                                                                               
  |======================================================================                                   |  66%
  |                                                                                                               
  |======================================================================                                   |  67%
  |                                                                                                               
  |=======================================================================                                  |  67%
  |                                                                                                               
  |=======================================================================                                  |  68%
  |                                                                                                               
  |========================================================================                                 |  69%
  |                                                                                                               
  |=========================================================================                                |  69%
  |                                                                                                               
  |=========================================================================                                |  70%
  |                                                                                                               
  |==========================================================================                               |  70%
  |                                                                                                               
  |==========================================================================                               |  71%
  |                                                                                                               
  |===========================================================================                              |  71%
  |                                                                                                               
  |===========================================================================                              |  72%
  |                                                                                                               
  |============================================================================                             |  72%
  |                                                                                                               
  |============================================================================                             |  73%
  |                                                                                                               
  |=============================================================================                            |  73%
  |                                                                                                               
  |=============================================================================                            |  74%
  |                                                                                                               
  |==============================================================================                           |  74%
  |                                                                                                               
  |==============================================================================                           |  75%
  |                                                                                                               
  |===============================================================================                          |  75%
  |                                                                                                               
  |================================================================================                         |  76%
  |                                                                                                               
  |================================================================================                         |  77%
  |                                                                                                               
  |=================================================================================                        |  77%
  |                                                                                                               
  |=================================================================================                        |  78%
  |                                                                                                               
  |==================================================================================                       |  78%
  |                                                                                                               
  |===================================================================================                      |  79%
  |                                                                                                               
  |====================================================================================                     |  80%
  |                                                                                                               
  |=====================================================================================                    |  81%
  |                                                                                                               
  |======================================================================================                   |  82%
  |                                                                                                               
  |=======================================================================================                  |  83%
  |                                                                                                               
  |========================================================================================                 |  83%
  |                                                                                                               
  |========================================================================================                 |  84%
  |                                                                                                               
  |=========================================================================================                |  84%
  |                                                                                                               
  |=========================================================================================                |  85%
  |                                                                                                               
  |==========================================================================================               |  85%
  |                                                                                                               
  |==========================================================================================               |  86%
  |                                                                                                               
  |===========================================================================================              |  86%
  |                                                                                                               
  |===========================================================================================              |  87%
  |                                                                                                               
  |============================================================================================             |  87%
  |                                                                                                               
  |============================================================================================             |  88%
  |                                                                                                               
  |=============================================================================================            |  88%
  |                                                                                                               
  |=============================================================================================            |  89%
  |                                                                                                               
  |==============================================================================================           |  89%
  |                                                                                                               
  |==============================================================================================           |  90%
  |                                                                                                               
  |===============================================================================================          |  90%
  |                                                                                                               
  |===============================================================================================          |  91%
  |                                                                                                               
  |================================================================================================         |  91%
  |                                                                                                               
  |================================================================================================         |  92%
  |                                                                                                               
  |=================================================================================================        |  92%
  |                                                                                                               
  |=================================================================================================        |  93%
  |                                                                                                               
  |==================================================================================================       |  93%
  |                                                                                                               
  |===================================================================================================      |  94%
  |                                                                                                               
  |===================================================================================================      |  95%
  |                                                                                                               
  |====================================================================================================     |  95%
  |                                                                                                               
  |====================================================================================================     |  96%
  |                                                                                                               
  |=====================================================================================================    |  96%
  |                                                                                                               
  |======================================================================================================   |  97%
  |                                                                                                               
  |=======================================================================================================  |  98%
  |                                                                                                               
  |======================================================================================================== |  99%
  |                                                                                                               
  |=========================================================================================================| 100%

snapshotDate(): 2021-05-18
see ?celldex and browseVignettes('celldex') for documentation
loading from cache
see ?celldex and browseVignettes('celldex') for documentation
loading from cache

We predict cell types by using SingleR function. It will predict cell types by using library from MonacoImmuneData

monaco.main

DataFrame with 10782 rows and 5 columns
                                                scores first.labels       tuning.scores          labels
                                              <matrix>  <character>         <DataFrame>     <character>
X112385.TCCGTCTTGATC    0.275612:0.176577:0.258834:...  Progenitors 0.1895553:0.1204320     Progenitors
X112385.AGTTTCCCTGGN    0.219137:0.134400:0.229646:...  Progenitors 0.0704632:0.0695319    CD4+ T cells
X112385.TGTTCCCATTCN    0.226374:0.158263:0.212045:...  Progenitors 0.0626434:0.0304610 Dendritic cells
X112385.CGGGACACTCCC    0.221196:0.180608:0.230647:...  Progenitors 0.0689471:0.0445722     Progenitors
X112385.ACCGGGTTACGG    0.223051:0.148000:0.223892:...  Progenitors 0.1701936:0.0990532     Progenitors
...                                                ...          ...                 ...             ...
X120386.TACAAGATATAC 0.1496382:0.0879767:0.1450494:...  Progenitors 0.0327718:0.0193404     Progenitors
X120386.CGCGATCAAATC 0.1689552:0.0836691:0.1564564:...  Progenitors 0.0327721:0.0249775         B cells
X120386.AGACTGGCTTTT 0.0821739:0.0865970:0.0898802:... CD4+ T cells 0.0456426:0.0356133       Basophils
X120386.CGATGGCCGAAT 0.1463504:0.0867209:0.1362387:...      B cells 0.0639621:0.0637729         B cells
X120386.AACGGTGCATAA 0.1531989:0.1175780:0.1478434:...  Progenitors 0.0885096:0.0770091     Progenitors
                       pruned.labels
                         <character>
X112385.TCCGTCTTGATC     Progenitors
X112385.AGTTTCCCTGGN    CD4+ T cells
X112385.TGTTCCCATTCN Dendritic cells
X112385.CGGGACACTCCC     Progenitors
X112385.ACCGGGTTACGG     Progenitors
...                              ...
X120386.TACAAGATATAC     Progenitors
X120386.CGCGATCAAATC         B cells
X120386.AGACTGGCTTTT       Basophils
X120386.CGATGGCCGAAT         B cells
X120386.AACGGTGCATAA     Progenitors

crc_dat_filter

class: SingleCellExperiment 
dim: 23768 10782 
metadata(0):
assays(4): counts log2counts cpm logcounts
rownames(23768): RP11-34P13.7 AL627309.1 ... KDM5D EIF1AY
rowData names(3): SYMBOL entrezid geneid
colnames(10782): X112385.TCCGTCTTGATC X112385.AGTTTCCCTGGN ... X120386.CGATGGCCGAAT
  X120386.AACGGTGCATAA
colData names(23): X batch_gid_dose ... label cell_type
reducedDimNames(2): PCA UMAP
mainExpName: NULL
altExpNames(0):

data.frame(colData(crc_dat_filter))

# colLabels(crc_dat_filter) <- factor(clust.kmeans$cluster)
plotReducedDim(crc_dat_filter, "UMAP", colour_by="cell_type")

Next, we would like to do DEseq analysis for cluster/cell_type analysis. Here, we are going to

• We need to create the count matrix and metadata for the counts

Create count matrix for each cluster/cell-type. Here, I am interested in the cell_types

nc

[1] 9

Getting number of cell per sample

table(crc_dat_filter$batch_gid_dose)

  112385_HCT116_0      112385_RKO_0    112385_SW480_0  112386_HCT116_10     112386_RKO_10   112386_SW480_10 
              294                83               382               132               115               187 
  112389_HCT116_0      112389_RKO_0    112389_SW480_0  112390_HCT116_10     112390_RKO_10   112390_SW480_10 
              216               272               276               367               517               526 
 120381_HCT116_50      120381_RKO_0  120381_SW480_200 120382_HCT116_200     120382_RKO_50    120382_SW480_0 
               25               695               474                28               291               657 
  120383_HCT116_0    120383_RKO_200   120383_SW480_50   120384_HCT116_0      120384_RKO_0    120384_SW480_0 
              721               142               627               780               623               648 
 120385_HCT116_50     120385_RKO_50   120385_SW480_50 120386_HCT116_200    120386_RKO_200  120386_SW480_200 
               92               328               690                15               117               462 

library(Matrix.utils)
# Subset metadata to only include the cluster and sample IDs to aggregate across
groups <- colData(crc_dat_filter)[, c("cell_type", "batch_gid_dose")]

pb <- aggregate.Matrix(t(counts(crc_dat_filter)), 
                       groupings = groups, fun = "sum") 

class(pb)

[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

dim(pb)

[1]   237 23768

pb[1:6,1:6]

6 x 6 sparse Matrix of class "dgCMatrix"
                         RP11-34P13.7 AL627309.1 RP11-34P13.13 AP006222.2 RP5-857K21.3 RP11-206L10.3
B cells_112385_HCT116_0             .          .             .          .            .             .
B cells_112385_RKO_0                .          .             .          .            .             .
B cells_112385_SW480_0              .          .             .          .            .             .
B cells_112386_HCT116_10            .          .             .          .            .             .
B cells_112386_RKO_10               .          .             .          .            .             .
B cells_112386_SW480_10             .          .             .          .            .             .

We will split our data by cell type; however, not always do all samples contain cells of every cell type. To determine which samples are present for each cell type we can run the following:

# Not every cluster is present in all samples; create a vector that represents how to split samples
splitf <- sapply(stringr::str_split(rownames(pb), 
                                    pattern = "_",  
                                    n = 2), 
                 `[`, 1)

# Turn into a list and split the list into components for each cluster and transform, so rows are genes and columns are samples and make rownames as the sample IDs
library(magrittr)
pb <- split.data.frame(pb, 
                       factor(splitf)) %>%
        lapply(function(u) 
                set_colnames(t(u), 
                             stringr::str_sub(rownames(u), str_locate(rownames(u),"_")[1]+1, -1)))

class(pb)

[1] "list"

# Explore the different components of list
str(pb)

List of 9
 $ B cells        :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:195960] 16 22 29 30 34 38 43 45 46 48 ...
  .. ..@ p       : int [1:30] 0 5639 11965 20625 25621 33172 40667 44163 49493 55525 ...
  .. ..@ Dim     : int [1:2] 23768 29
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:29] "112385_HCT116_0" "112385_RKO_0" "112385_SW480_0" "112386_HCT116_10" ...
  .. ..@ x       : num [1:195960] 4 14 1 1 2 2 1 15 1 17 ...
  .. ..@ factors : list()
 $ Basophils      :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:13711] 16 48 111 126 149 158 159 166 181 231 ...
  .. ..@ p       : int [1:12] 0 1208 1772 3101 4002 5559 6534 6861 9007 11444 ...
  .. ..@ Dim     : int [1:2] 23768 11
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:11] "112386_HCT116_10" "112389_HCT116_0" "112390_HCT116_10" "120381_HCT116_50" ...
  .. ..@ x       : num [1:13711] 1 1 1 2 1 1 1 2 1 1 ...
  .. ..@ factors : list()
 $ CD4+ T cells   :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:190457] 16 21 22 25 30 34 39 43 45 46 ...
  .. ..@ p       : int [1:31] 0 6919 13425 19828 24505 32357 36334 41247 48810 53854 ...
  .. ..@ Dim     : int [1:2] 23768 30
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:30] "112385_HCT116_0" "112385_RKO_0" "112385_SW480_0" "112386_HCT116_10" ...
  .. ..@ x       : num [1:190457] 11 2 15 1 1 2 2 1 21 5 ...
  .. ..@ factors : list()
 $ CD8+ T cells   :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:97781] 39 45 48 73 126 158 159 166 184 197 ...
  .. ..@ p       : int [1:28] 0 1546 2829 6699 8811 12793 15052 17306 20677 23896 ...
  .. ..@ Dim     : int [1:2] 23768 27
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:27] "112385_HCT116_0" "112385_RKO_0" "112385_SW480_0" "112386_HCT116_10" ...
  .. ..@ x       : num [1:97781] 1 2 1 1 1 2 1 5 1 1 ...
  .. ..@ factors : list()
 $ Dendritic cells:Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:241222] 10 16 17 21 22 26 29 30 31 34 ...
  .. ..@ p       : int [1:31] 0 9543 17305 27212 33558 39695 47993 51492 56234 62383 ...
  .. ..@ Dim     : int [1:2] 23768 30
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:30] "112385_HCT116_0" "112385_RKO_0" "112385_SW480_0" "112386_HCT116_10" ...
  .. ..@ x       : num [1:241222] 1 28 1 5 93 1 7 16 1 7 ...
  .. ..@ factors : list()
 $ Monocytes      :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:97970] 30 45 47 55 60 91 111 121 134 135 ...
  .. ..@ p       : int [1:30] 0 1752 4311 6119 8356 11396 12532 14893 16932 18244 ...
  .. ..@ Dim     : int [1:2] 23768 29
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:29] "112385_HCT116_0" "112385_RKO_0" "112385_SW480_0" "112386_HCT116_10" ...
  .. ..@ x       : num [1:97970] 2 1 1 1 1 2 1 1 1 1 ...
  .. ..@ factors : list()
 $ NK cells       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:42278] 66 166 193 197 315 353 365 412 450 469 ...
  .. ..@ p       : int [1:23] 0 560 3331 5095 6504 8356 9829 10469 11364 12979 ...
  .. ..@ Dim     : int [1:2] 23768 22
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:22] "112385_HCT116_0" "112385_SW480_0" "112386_SW480_10" "112389_HCT116_0" ...
  .. ..@ x       : num [1:42278] 2 2 1 1 1 1 5 2 1 2 ...
  .. ..@ factors : list()
 $ Progenitors    :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:318652] 1 11 15 16 18 20 21 22 23 25 ...
  .. ..@ p       : int [1:31] 0 13307 20844 34388 44654 51492 64013 72542 79845 89165 ...
  .. ..@ Dim     : int [1:2] 23768 30
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:30] "112385_HCT116_0" "112385_RKO_0" "112385_SW480_0" "112386_HCT116_10" ...
  .. ..@ x       : num [1:318652] 1 12 1 153 1 1 17 354 3 3 ...
  .. ..@ factors : list()
 $ T cells        :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. ..@ i       : int [1:110444] 11 16 38 45 48 91 116 124 126 132 ...
  .. ..@ p       : int [1:30] 0 2330 5071 8090 10208 11374 14027 17744 19437 22991 ...
  .. ..@ Dim     : int [1:2] 23768 29
  .. ..@ Dimnames:List of 2
  .. .. ..$ : chr [1:23768] "RP11-34P13.7" "AL627309.1" "RP11-34P13.13" "AP006222.2" ...
  .. .. ..$ : chr [1:29] "112385_HCT116_0" "112385_SW480_0" "112386_HCT116_10" "112386_RKO_10" ...
  .. ..@ x       : num [1:110444] 1 1 1 2 5 1 1 1 2 1 ...
  .. ..@ factors : list()

# Print out the table of cells in each cluster-sample group
options(width = 100)
table(crc_dat_filter$cell_type, crc_dat_filter$batch_gid_dose)

                 
                  112385_HCT116_0 112385_RKO_0 112385_SW480_0 112386_HCT116_10 112386_RKO_10
  B cells                      13           15             47               12            24
  Basophils                     0            0              0                2             0
  CD4+ T cells                 26           21             27               11            39
  CD8+ T cells                  2            1              7                3             7
  Dendritic cells              47           20             45               20            14
  Monocytes                     3            2              2                2             5
  NK cells                      1            0              3                0             0
  Progenitors                 199           24            245               78            22
  T cells                       3            0              6                4             4
                 
                  112386_SW480_10 112389_HCT116_0 112389_RKO_0 112389_SW480_0 112390_HCT116_10
  B cells                      29              17           35             42               28
  Basophils                     0               1            0              0                3
  CD4+ T cells                 12              25           94             26               24
  CD8+ T cells                  3               6           15             11                2
  Dendritic cells              17              16           21             42               54
  Monocytes                     2               7            5              3               13
  NK cells                      2               3            4              3                1
  Progenitors                 120             134           81            145              232
  T cells                       2               7           17              4               10
                 
                  112390_RKO_10 112390_SW480_10 120381_HCT116_50 120381_RKO_0 120381_SW480_200
  B cells                    82              50                0           62               49
  Basophils                   0               0                1            0                2
  CD4+ T cells              122              32                2          134               27
  CD8+ T cells               18              13                0           53               10
  Dendritic cells           102              63                5          124               69
  Monocytes                  25               5                2            7                3
  NK cells                    2               3                0            8                3
  Progenitors               143             348               13          265              297
  T cells                    23              12                2           42               14
                 
                  120382_HCT116_200 120382_RKO_50 120382_SW480_0 120383_HCT116_0 120383_RKO_200
  B cells                         2            50             64              16             23
  Basophils                       0             0              0               1              1
  CD4+ T cells                    1            61             42              36             22
  CD8+ T cells                    0            17             21              11              5
  Dendritic cells                10            76             66              92             30
  Monocytes                       0            24              9              30             15
  NK cells                        0             2              5               1              1
  Progenitors                    14            44            418             507             39
  T cells                         1            17             32              27              6
                 
                  120383_SW480_50 120384_HCT116_0 120384_RKO_0 120384_SW480_0 120385_HCT116_50
  B cells                      33              35           78             81                1
  Basophils                     7               5            0              0                0
  CD4+ T cells                 23              57          164             48                7
  CD8+ T cells                  9               5           15             14                1
  Dendritic cells              33             101           89             87               23
  Monocytes                     7              16           14              7                3
  NK cells                      9               3            4              7                0
  Progenitors                 487             545          243            390               52
  T cells                      19              13           16             14                5
                 
                  120385_RKO_50 120385_SW480_50 120386_HCT116_200 120386_RKO_200 120386_SW480_200
  B cells                    59              41                 1             22               57
  Basophils                   1               0                 0              0                3
  CD4+ T cells               45              24                 1             16               28
  CD8+ T cells                7               7                 0              6                8
  Dendritic cells           110              85                 2             33               78
  Monocytes                  24               4                 1              7                8
  NK cells                    1               5                 0              0                1
  Progenitors                66             518                 9             27              268
  T cells                    15               6                 1              6               11

Sample-level metadata

# Get sample names for each of the cell type clusters

# prep. data.frame for plotting
get_sample_ids <- function(x){
        pb[[x]] %>%
                colnames()
}

de_samples <- map(1:nc, get_sample_ids) %>%
        unlist()

pb[[3]]  %>%
                colnames()

 [1] "112385_HCT116_0"   "112385_RKO_0"      "112385_SW480_0"    "112386_HCT116_10"  "112386_RKO_10"    
 [6] "112386_SW480_10"   "112389_HCT116_0"   "112389_RKO_0"      "112389_SW480_0"    "112390_HCT116_10" 
[11] "112390_RKO_10"     "112390_SW480_10"   "120381_HCT116_50"  "120381_RKO_0"      "120381_SW480_200" 
[16] "120382_HCT116_200" "120382_RKO_50"     "120382_SW480_0"    "120383_HCT116_0"   "120383_RKO_200"   
[21] "120383_SW480_50"   "120384_HCT116_0"   "120384_RKO_0"      "120384_SW480_0"    "120385_HCT116_50" 
[26] "120385_RKO_50"     "120385_SW480_50"   "120386_HCT116_200" "120386_RKO_200"    "120386_SW480_200" 

# Get cluster IDs for each of the samples

samples_list <- purrr::map(1:nc, get_sample_ids)

get_cluster_ids <- function(x){
        rep(names(pb)[x], 
            each = length(samples_list[[x]]))
}

de_cluster_ids <- purrr::map(1:nc, get_cluster_ids) %>%
        unlist()


get_sample_ids(3)

 [1] "112385_HCT116_0"   "112385_RKO_0"      "112385_SW480_0"    "112386_HCT116_10"  "112386_RKO_10"    
 [6] "112386_SW480_10"   "112389_HCT116_0"   "112389_RKO_0"      "112389_SW480_0"    "112390_HCT116_10" 
[11] "112390_RKO_10"     "112390_SW480_10"   "120381_HCT116_50"  "120381_RKO_0"      "120381_SW480_200" 
[16] "120382_HCT116_200" "120382_RKO_50"     "120382_SW480_0"    "120383_HCT116_0"   "120383_RKO_200"   
[21] "120383_SW480_50"   "120384_HCT116_0"   "120384_RKO_0"      "120384_SW480_0"    "120385_HCT116_50" 
[26] "120385_RKO_50"     "120385_SW480_50"   "120386_HCT116_200" "120386_RKO_200"    "120386_SW480_200" 

get_cluster_ids(2)

 [1] "Basophils" "Basophils" "Basophils" "Basophils" "Basophils" "Basophils" "Basophils" "Basophils" "Basophils"
[10] "Basophils" "Basophils"

names(pb)[2]

[1] "Basophils"

Subsetting dataset to cell(s) of interest

# Generate vector of cluster IDs
clusters <- levels(factor(metadata$cell_type))
clusters

[1] "B cells"         "Basophils"       "CD4+ T cells"    "CD8+ T cells"    "Dendritic cells" "Monocytes"      
[7] "NK cells"        "Progenitors"     "T cells"        

# Subset the metadata to only the B cells
cluster_metadata <- metadata[which(metadata$cell_type == clusters[1]), ]
head(cluster_metadata)

# Assign the rownames of the metadata to be the sample IDs
rownames(cluster_metadata) <- cluster_metadata$batch_gid_dose
head(cluster_metadata)

# Run DESeq2 differential expression analysis
dds <- DESeq(dds)

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
-- replacing outliers and refitting for 1 genes
-- DESeq argument 'minReplicatesForReplace' = 7 
-- original counts are preserved in counts(dds)
estimating dispersions
fitting model and testing

# Plot dispersion estimates
plotDispEsts(dds)

counts[, which(colnames(counts) %in% rownames(cluster_metadata))]

23768 x 29 sparse Matrix of class "dgCMatrix"

   [[ suppressing 29 column names ‘112385_HCT116_0’, ‘112385_RKO_0’, ‘112385_SW480_0’ ... ]]
   [[ suppressing 29 column names ‘112385_HCT116_0’, ‘112385_RKO_0’, ‘112385_SW480_0’ ... ]]

                                                                               
RP11-34P13.7   . .  . . . . . . . . . . .  . . .  . . .  .  .  . . . .  . . . .
AL627309.1     . .  . . . . . . . . . . .  . . .  . . .  .  .  . . . .  . . . .
RP11-34P13.13  . .  . . . . . . . . . . .  . . .  . . .  .  .  . . . .  . . . .
AP006222.2     . .  . . . . . . . . . . .  . . .  . . .  .  .  . . . .  . . . .
RP5-857K21.3   . .  . . . . . . . . . . .  . . .  . . .  .  .  . . . .  . . . .
RP11-206L10.3  . .  . . . . . . . . . . .  . . .  . . .  .  .  . . . .  . . . .
RP11-206L10.2  . .  . . . . . . . 1 . 2 .  . . .  . . .  .  .  . . . .  . . . .
RP11-206L10.9  . .  . . . . . . . . . . .  3 . .  . . .  .  .  . . . 1  1 . . 1
RP11-206L10.10 . .  . . . . . . . . . 1 .  . . .  . . .  .  .  . . . .  . . . .
FAM87B         . .  . . . . . . . . . . .  . . .  . . .  .  .  . . . .  . . . .
LINC00115      . .  . . . 1 . . . 1 . . .  . . .  1 . 1  .  .  . . . .  . . . .
LINC01128      . .  3 1 . . 1 . . . 3 1 2  1 . 1  1 . .  .  1  . 1 . .  . . . .
RP11-54O7.16   . .  . . . . . . . . . . .  . . .  . . .  .  .  . . . .  . . . .
RP11-54O7.1    . .  . . . . . . . . . . .  . . .  . . .  .  .  . 1 . .  . . . .
RP11-54O7.3    . .  1 . . . . . . . . . .  . . .  1 . .  .  .  . . . .  . . . .
SAMD11         . 1  2 . . 3 . . . . . 3 .  . . .  9 . 1 11  .  . 5 . .  6 . . .
NOC2L          4 5 12 1 2 2 2 1 2 4 7 5 9 13 . 2 16 2 4  8 11 11 9 . 7 11 . 2 7

 ..............................
 ........suppressing 23734 rows in show(); maybe adjust 'options(max.print= *, width = *)'
 ..............................

   [[ suppressing 29 column names ‘112385_HCT116_0’, ‘112385_RKO_0’, ‘112385_SW480_0’ ... ]]

                                                                        
CLIC2         . . . . . . . . . . . . .  . . . . . . . . . . . . . . . .
RP13-228J13.1 . . . . . . . . . . . . .  . . . . . . . . 1 . . . . . . .
F8A3          . . . . . . . . . . . . .  . . . . . . . . . . . . . . . .
TMLHE-AS1     . . . . 1 . . . . . . . .  . . . . . . . . . . . . . . . .
RP11-218L14.4 . . . . . . . . . . . . .  . . . . . . . . . . . . . . . .
TMLHE         . . 8 . 4 4 1 . 2 1 2 1 2  3 . 1 . . 1 . . . 3 . 5 2 . . 7
SPRY3         . . . . . . . . . . . . .  . . . . . . . . . . . . . . . .
VAMP7         . 2 . 1 1 5 . . . 2 . 3 1 11 . 1 2 . 2 2 1 1 3 . 3 3 . 4 8
AJ271736.10   . . . . . . . . . 1 . . .  1 . . . . . . . . . . . 3 . . .
WASIR1        . . . . . . . . . . . . .  . . . . . . . . 1 . . . . . . .
RPS4Y1        . . . . . . . . . . . . .  . . . . . . . . . . . . . . . .
AC012667.1    1 . . . . 1 . . . . . . .  . . . . . . . 2 1 . . . . . . .
TMSB4Y        . . . . . . . . . . . . .  . . . . . . . . . . . . . . . .
NLGN4Y        . . . . . . . . . . . . .  . . . . . . . . . . . . 1 . . .
TTTY14        . . . . . . . . . . . . .  . . . . . . . . . . . . . . . .
KDM5D         . . . . . . . . . . . . .  . . . . . . . . . . . . . . . .
EIF1AY        . . . . . . . 1 . . 1 . .  . . 1 . . 1 . . . . . 1 . . . .

Results from DESeq

resultsNames(dds)

[1] "Intercept"              "df.gid_RKO_vs_HCT116"   "df.gid_SW480_vs_HCT116"

res <- lfcShrink(dds, 
                 coef = 2,
                 res=res)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895

# Turn the results object into a tibble for use with tidyverse functions
res_tbl <- res %>%
        data.frame() %>%
        rownames_to_column(var="gene") %>%
        as_tibble()

# Check results output
res_tbl

# Write all results to file
# write.csv(res_tbl,
#           paste0("results/", clusters[1], "_", levels(cluster_metadata$sample)[2], "_vs_", levels(cluster_metadata$sample)[1], "_all_genes.csv"),
#           quote = FALSE, 
#           row.names = FALSE)

Look for significant genes

# Set thresholds
padj_cutoff <- 0.05

# Subset the significant results
sig_res <- dplyr::filter(res_tbl, padj < padj_cutoff) %>%
        dplyr::arrange(padj)

# Check significant genes output
sig_res